// MLE for predicted Wt. proportion in the diet of pred j
DATA_SECTION
  // read in data from predWt.dat
  init_int LWnobs
  init_int Nnobs
  init_int wtdatnobs
  init_int maxdatnobs
  init_int npars

  init_vector Lp(1,LWnobs)
  init_vector Wp(1,LWnobs)
  init_vector prey_size(1,Nnobs)
  init_vector Np(1,Nnobs)
  init_vector propN(1,Nnobs)
  init_vector Wt(1,Nnobs)
  init_vector pred_size(1,wtdatnobs)
  init_vector wt(1,wtdatnobs)
  init_vector G(1,wtdatnobs)
 LOCAL_CALCS
   for (int l=1;l<=wtdatnobs;l++)
   {
     if( G(l)==0)
     {
       G(l)=0.0001;
     }
     if( G(l)==1)
     {
       G(l)=0.9999;
     }
   }
 END_CALCS
  init_vector WtTot(1,wtdatnobs)
  init_vector juse(1,maxdatnobs)
  init_vector z(1,maxdatnobs)
  init_vector pars(1,npars)
  !! cout << pars<<endl;
  
  int l // pre-allocate k
  int i // ditto
  int j //ditto
  
INITIALIZATION_SECTION
  alphap 1.89
  betap 0.09
  sigmaw 3.2
  beta0g -2.6
  beta1g 1.86
  beta2g 0.24
  beta3g 0.0018
  beta4g -0.66
  logphig -0.9
  q 3.4
  alphaz 1.9
  betaz .17
  logshapez -1.7

PARAMETER_SECTION
  number bb
  vector aa(1,maxdatnobs)
  number shapez // pre-allocate k
  number phig   // ditto
  init_number logshapez // pre-allocate k
  init_number logphig // ditto
  
  init_number alphap
  init_number betap
  init_number sigmaw
  init_number beta0g
  init_number beta1g
  init_number beta2g
  init_number beta3g
  init_number beta4g
  init_number q
  init_number alphaz
  init_number betaz
  
  vector logWphat(1,LWnobs)
  vector Wphat(1,LWnobs)
  vector epsw(1,LWnobs)
  vector f(1,LWnobs)
  vector NLLw(1,LWnobs)
  vector logzmean(1,maxdatnobs)
  vector zmean(1,maxdatnobs)
  vector NLLz(1,maxdatnobs)
  vector limhat(1,maxdatnobs)
  vector lambda(1,wtdatnobs)
  vector logitphi(1,Nnobs)
  vector phi(1,wtdatnobs)
  vector mu(1,wtdatnobs)
  vector logitmu(1,wtdatnobs)
  vector a(1,wtdatnobs)
  vector b(1,wtdatnobs)
  vector Guse(1,wtdatnobs)
  vector NLLg(1,wtdatnobs)
  objective_function_value obj_fun;

//PRELIMINARY_CALCS_SECTION
  
PROCEDURE_SECTION
  shapez = mfexp(logshapez);
  phig   = mfexp(logphig);
//  calc_LW();
//  cal_lim();
//  cal_lambda();
//  calc_wtprop();

//FUNCTION calc_LW(alphap betap)
  // calculate the lw weight relatioships
  logWphat = alphap+Lp*betap;
  epsw     = log(Wp)-logWphat;
  obj_fun +=-0.5*norm2(epsw); // double check that this is correct - need to be log normal errors....
  //f=(norm2(epsw)); 
 // obj_fun +=LWnobs/2.*log(f);    // make it a likelihood function?

//FUNCTION calc_lim(alphaz betaz shapez)
  //calulate the max size prey available for each pred size
  logzmean = alphaz+betaz*juse;
  zmean    = mfexp(logzmean);
  bb       = shapez;  // scalar
  aa       = shapez*zmean; // vector
  //for (i=1;i<=maxdatnobs;i++)
  // neg log lik for gamma dist
  NLLz     = -1*(aa*log(bb)-gammln(aa)+elem_prod((aa-1),log(zmean))+-bb*zmean);
  obj_fun +=sum(NLLz);
  
//FUNCTION calc_lambda(alphaz betaz alphap betap q)
  for (j=1;j<=Nnobs;j++)
  {
    limhat(j) = alphaz+betaz*pred_size(j); // scalar based on size of predator (j)
    // Wphat     = alphap*pow(prey_size,betap);  // vector of predicted prey weights
    // logitphi  = elem_div(limhat(j)*q-prey_size,prey_size); // vector of phis
    logitphi(j)  = elem_div(limhat(j)*q-prey_size(j),prey_size(j)); // vector of phis
    phi       = elem_div(mfexp(logitphi),1+mfexp(logitphi));
    for (i=1;i<=maxdatnobs;i++) 
    {
      if ( prey_size(i)==0)
      {
         phi(i)=1;
       }
    }
    //lambda(j)=sum(elem_prod(elem_prod(Np,Wphat),phi))/sum(elem_prod(Np,Wphat));    // scalar based on biomass
    lambda(j) = sum(elem_prod(Np,phi))/sum(Np); //          # numerically based
  }
  
//FUNCTION calc_wtprop(beta0g beta1g beta2g beta3g beta4g phig)
  logitmu = beta0g+beta1g*lambda+beta2g*pred_size+beta3g*elem_prod(pred_size,pred_size)+beta4g*elem_prod(pred_size,lambda);
  mu      = elem_div(mfexp(mu),1+mfexp(mu));
  a       = mu*phig; // vector
  b       = (1.-mu)*phig;  // vector
  // neg log lik for beta dist
  NLLg=-1*(gammln(a+b)-gammln(a)+gammln(b)+elem_prod(a-1,log(mu))+elem_prod(b-1,log(1-mu)));
  // NLL.g[is.na(NLL.g)]=1000000 penalize function?
  obj_fun +=sum(NLLg);
  
  //sum(sum.na(NLL.z),sum.na(NLL.w),sum.na(NLL.g))
 
REPORT_SECTION
  cout <<endl<<"========End of phase: "<<current_phase()<<" ============"<<endl<<endl;

  report << "predicted proportion of wt in the stomach"<<endl;
  report << mu <<endl;



TOP_OF_MAIN_SECTION
  arrmblsize = 10000000; 
